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Abstract 

Ultra-high energy cosmic rays generate extensive air showers in Earth's atmosphere. 
A standard approach to reconstruct the energy of an ultra-high energy cosmic rays 
is to sample the lateral profile of the particle density on the ground of the air shower 
with an array of surface detectors. 

For cosmic rays with large inclinations, this reconstruction is based on a model 
of the lateral profile of the muon density observed on the ground, which is fitted 
to the observed muon densities in individual surface detectors. The best models 
for this task are derived from detailed Monte-Carlo simulations of the air shower 
development. We present a phenomenological parametrization scheme which allows 
to derive a model of the average lateral profile of the muon density directly from a 
fit to a set of individual Monte-Carlo simulated air showers. The model reproduces 
the detailed simulations with a high precision. As an example, we generate a muon 
density model which is valid in the energy range 10^^ eV < < 10^0 eV and the 
zenith angle range 60° < 9 < 90°. 

We will further demonstrate a way to speed up the simulation of such muon 
profiles by three orders of magnitude, if only the muons in the shower are of interest. 
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1 Introduction 



Ultra-high energy cosmic rays (UHECRs) are cosmic rays with energies above 
1 EeV = 10^*^ eV. They have been under study for several decades, still their 
origin is not well known. Information about the origin is encoded in the energy 
spectrum [1] , the mass composition [2] and a possible anisotropy of their arrival 
directions |3], which can be measured experimentally. 

The huge energy and the low flux (roughly 1 particle per km^ per year above 10 
EeV) make a direct measurement of the momentum and energy of a UHECR 
with balloon or satellite experiments unfeasible. Instead, they are observed 
indirectly with ground-based detectors that use Earth's atmosphere as a large 
calorimeter. The interactions of the UHECR with atmospheric nuclei generate 
an extensive particle shower which is sampled by these detectors. 

One realization of such a detector is a large surface array of particle counters. 
The array samples the lateral profile of the particle density of the shower which 
consist mainly of photons, electrons, and muons. The arrival direction of the 
UHECR can be obtained rather directly from the measured arrival time of the 
shower front in individual particle counters. The reconstruction of the energy 
E of the cosmic ray is more complex and requires to fit a model of the lateral 
particle distribution around the shower axis to the measured particle counts. 

Most surface array experiments concentrate on showers with inclinations less 
than 60°. For such showers, the particle distribution is radially symmetric in 
good approximation and well described by comparably simple empirical mod- 
els of the NKG-type [HIS]. At larger inclinations, the effect of the geomagnetic 
field on the particle distribution cannot be neglected. The symmetry becomes 
broken and the NKG-type models fail to describe the particle distribution. 

In an ideal surface array, these very inclined showers constitute 25 % of the 
number of arriving events. Recovering them yields a significant gain in the 
event statistics of the experiment. The ability to reconstruct very inclined 
showers also increases the field of view of the detector and thus the total ob- 
servable region of the sky. This is particularly relevant for anisotropy searches. 

It was first demonstrated with the Water- Cherenkov detectors of the Haverah 
Park experiment [B] that the energy E of cosmic rays at large zenith angles 
6 > 60° can be derived from the total number of muons N^j^ which arrive 
at the ground [3 [8]. The same approach is now used by the Pierre Auger 
Observatory [9]. The number of muons A*"^ on the ground is obtained from a 
fit of a model of the average lateral profile of the muon densit}(l] on the 



The term "muon density" in this article refers to the time-integrated particle flux 
through the ground plane initiated by an air shower. 
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ground to the measured signals. This fit exploits the following factorization 



n^~iV^(E,A^)x^(x,i/,^,0), (1) 

whereas Nfj_ is the number of muons on the ground which depends only on 
the energy E, mass A, and inclination 6 of the cosmic ray, while is a 
normalized lateral profile of the muon density which depends only on the 
ground coordinates {x,y) and the shower direction {0,(j)). The normalized 
profile also depends on the properties of the observation site like the ground 
altitude, the geomagnetic field and the atmosphere, but those are considered 
fixed here. In simple terms, the factorization says that the shape of the lateral 
profile remains the same for all showers arriving from a certain direction in 
very good approximation, while its amplitude carries all the information about 
the energy E and mass A of the cosmic ray. This invariance of the profile shape 
is called shower universality. 

The universality is very useful, because is too complex to be fitted from 
the data sampled by the surface array on an event-by-event basis. It has to 
be modeled. However, if is predicted by a model, the reconstruction of 
reduces to a fit of three parameters: the two intersection coordinates of the 
shower axis with the ground and the amplitude Nfj_. 

This reconstruction approach even works if the particle counters of the sur- 
face detector cannot distinguish between different species of charged particles. 
Very inclined air showers arrive in a very late stage of their development on 
the ground, where the only remaining electromagnetic particles in the shower 
are generated by the muons themselves, mostly via decay. Therefore, signals 
generated by such old showers remain proportional to the local muon den- 
sity n^, because the electromagnetic particles only enhance the signal by a 
constant factor in first approximation [TU] . 

Models of the normalized profile are currently derived from detailed Monte- 
Carlo simulations of extensive air showers, which represent best our current 
theoretical knowledge. However, it is not feasible to make a Monte- Carlo sim- 
ulation of for every possible shower direction (6',0), simply because these 
simulations consume considerable computing time and and storage space. The 
established solution to this issue is a semi-analytical model of which is 
based on detailed simulations but predicts the azimuthal dependency of 
analytically |7]. The semi-analytical model is instructive and reproduces the 
detailed simulations well, but still deviates somewhat from full air shower sim- 
ulations since some effects are neglected in the analytical part of the model. 

In this article, we follow a different path and propose a phenomenological 
parametrization of which can be fitted directly to a set of simulated air 
showers. The parametrization is not derived from an analytical theory, it only 
exploits some general principles of the air shower development. The approach 
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is valid up to at least 4 km from the shower axis in the energy range 10^^ eV < 
E < 10^° eV and the zenith angle range 60° < 9 < 88° and reproduces the 
output of detailed simulations better than the semi-analytical model. As a 
consequence, our model should lead to smaller biases in the reconstructed 
muon number on the ground. 

The parametrization procedure does not require a specific distribution of sim- 
ulated showers over the parameter range {9, cj)) of interest as input as long as 
that range is sufficiently covered with simulated showers. It can therefore be 
applied to many existing air shower libraries. As a consequence of the shower 
universality, the phenomenological parametrization of is practically inde- 
pendent of the nature of the primary particle, except for a global factor. 

In this article, we derive and motivate the parametrization. It is then applied to 
a set of 1800 simulated air showers as an example and in order to demonstrate 
its precision. On a related side note, we point out an efficient way to speed 
up the detailed simulation of n^-profiles, which can then be used as input 
of the parametrization. In particular, we show that exploiting the azimuthal 
symmetry of the muon profile at the production point and neglecting the 
calculation of the electromagnetic cascade leads to computation times much 
smaller than conventional Monte-Carlo simulations. 

The article is organized as follows. In Section |2l we review some general fea- 
tures of very inclined air showers in order to motivate the parametrization 
approach. The library of air showers used for the example application and the 
fast simulation approach are described in Section 121 The parametrization is 
then discussed in Section H] and its precision demonstrated. A summary of the 
results is given in Section |5l 



2 Very inclined air showers: the general picture 

The physics (see e.g. ref. [TtlSI ITTHTB] ) which influence the muon component of 
very inclined air showers are reviewed in this section. 

The lateral proflle of the muon density on the ground has different prop- 
erties for standard showers (0° < 6' < 60°) and for very inclined showers 
(60° < 9 < 90°). Standard showers are dominated by photons and electrons 
generated in the hadronic interactions. The flux of particles through the shower 
front plane is almost radially symmetric with respect to the shower axis. 

Very inclined showers travel longer through the atmosphere. As the shower 
arrives at the ground, the electromagnetic component generated by hadronic 
interactions is already absorbed. What remains is a muon shower in first ap- 
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proximation. The lateral profile of the muon density gets deformed by the 
geomagnetic field, atmospheric attenuation, and geometrical effects, so that 
the particle flux through the shower front plane is no longer radially symmet- 
ric. These asymmetries become very large as 9 approaches 90°. 

The primary electromagnetic component coming from the decay of neutral pi- 
ous is already almost absorbed when the shower reaches the ground a,i 9 60° 
and can be considered extinct at 6* > 70° (see also Fig. [3]). A non-negligible 
electromagnetic component is still detectable on the ground at all angles which 
is produced by the muons themselves. The muons produce photons and elec- 
trons mostly through decay, but also via bremsstrahlung, e"'"e~-production and 
delta rays along their path through the atmosphere. In some sense, they are 
surrounded by electromagnetic sub-showers with a moderate lateral extension 
- typically a few tens of meters. 

Therefore, as long as the muon density is larger than about 0.1 m~^, these 
electromagnetic sub-showers overlap result in a continuous halo. The elec- 
tromagnetic particles arrive together with the muons within a few tens of 
nanoseconds and do not increase the longitudinal thickness of the shower 
front. However, they have to be considered in particle detectors that can- 
not distinguish between muons and the electromagnetic halo. We will ignore 
the impact of electromagnetic particles in the following. We assume that ei- 
ther the detectors distinguish between different particle species, so that the 
electromagnetic component can be ignored, or that the close proportionality 
between electromagnetic and muonic components in very inclined showers is 
exploited in the data analysis in order to estimate the muon density from the 
mixed signal [5|[TU]. 



2.1 Conventions 



Quantitative calculations and detailed simulations of extensive air showers 
depend on the observation site. Important are the local geomagnetic field 
the altitude of the ground hg above sea level and the air density profile of 
the air above the site. In this article, we do all calculations for the site of the 
Southern Pierre Auger Observatory in Malargiie, Argentina. 

The ground plane altitude is taken as 1425 m. The geomagnetic field B is 
treated as a constant fielcOJ: 

B = 24.6 /iT, 5b = 4.2°, 9b = -35.2°, 



^ The geomagnetic field currently varies by about 1° in direction and 2 % in mag- 
nitude over 10 years in Malargiie [17j . 
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Fig. 1. The shower front plane coordinate system [7]: is anti-parallel to the mo- 
mentum vector of the shower, is parallel to Bt, the component of the geomagnetic 
field projected into the shower plane, r and ip are polar coordinates in the shower 
front plane. 

whereas 6b and 6b are the geomagnetic declination and inclination. 

An average profile of the air density over the site is approximated by the US 
standard atmosphere [18] . 



The muon profiles derived in this article are supposed to be comparable with 
experimental measurements and therefore need to regard the energy thresh- 
old of the applied particle counters. Scintillators detect practically all muons; 
their threshold for Cherenkov emission in water is about 50 MeV. However 
the threshold for buried or shielded detectors can reach a few GeV. The sim- 
ulations shown in this article are done for a muon energy threshold 

Efhres = 0.145 GeV. 

The results of this article remain qualitatively correct if it is increased up to a 
few GeV, but quantitatively they depend on the threshold. The dependency 
of muon profiles on the observation site and the muon energy threshold is not 
mentioned explicitly in the rest of the article. 

In order to discuss the lateral profile of an extensive air shower, it is useful 
to introduce a special coordinate system: the shower front plane coordinate 
system (see Fig.[T]). The shower front plane is perpendicular to the shower 
axis. All observations are still done in the ground plane, but the coordinates 
are projected onto this plane. The projection restores some of the principal 
symmetries of the shower profile. We will refer to this coordinate system during 
the rest of the article. 

Finally, we use the convention of the Pierre Auger Observatory for the az- 
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Fig. 2. As a function of the zenith angle at the ground altitude of 1425 m, we 
show: a) the integrated atmospheric slant depth X-^tra along the shower path, b) the 
distance d between the shower maximum and the impact of the core on the ground, 
and c) the altitude h of the electromagnetic shower maximum above the ground 
level. The latter are shown for different energies, because the depth X^ax of the 
shower maximum depends on the energy E. For the calculations of the slant depth 
X along a shower path in a curved atmosphere see e.g. [TT]. A parametrization of 
XinaxiE) is taken from [2]. 

imuthal angle 0, where a shower arriving from the geographic East equals to 
(/) = 0°, and a shower from the geographic North equals to = 90°. 



2.2 Development of the muon component 



An air shower induced by a proton or a nucleus first produces a hadronic 
cascade which mostly produces charged and neutral pions in each step. The 
neutral pions decay almost immediately into two photons and feed an electro- 
magnetic component. The decay of the charged pions feeds a muon component. 
The total atmospheric slant depth Xatm increases with the zenith angle 6; for 
example, for an observer at 1400 m above sea level, when 6 rises from 0° to 
90°, the depth down to the ground increases from 870gcm~^ to 31000 gcm~^ 
(see Fig.[2D. 

At an inclination 6 > 60°, the total atmospheric depth Xatm is more than twice 
the depth X^ax of the air shower maximum, which ranges between 650 g cm~^ 
and 800 gcm~^ between 10^^ eV and 10^*^ eV [2], depending on the energy E 
and mass A of the cosmic ray. Xatm is more than three times the depth of the 
end of the hadronic cascade, where most muons are produced (see Fig.[3]). As a 
consequence, the electromagnetic cascade is fully extinguished at ground level. 
Only the most energetic muons survive, accompanied by an electromagnetic 
halo produced by the decay and other radiative processes of the muons. 
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Fig. 3. Comparison of two inclined proton showers, simulated with CORSIKA using QGSJet-II and FLUKA for the hadronic inter- 
actions [22fl24j ■ The energies, zenith and azimuth angles of two showers are given in the legend right of the plots. The plots show a) 
the total particle number vs. vertical atmospheric depth, b) the particle number per tank on the ground vs. radial distance from the 
shower axis, and c) the momentum spectrum on the ground. In b), r > relates to particles collected within ±10° around the direction 
perpendicular to the projected geomagnetic field Bt and r < (to be taken in absolute value) relates to particles collected within 
±10° around the direction parallel to Bt- In all plots, only particles above certain energy thresholds are shown: -Ethres = 250 keV for 
electromagnetic particles and -Ethres = 0.1 GeV for everything else. The low energy photon peak in c) is caused by /e~ annihilation. 
The particle depletion along Bt in b) for 6 = 82° is the result of geomagnetic deflection. The rapid decay of the muon number at large 
slant depths is artificial and caused by the lateral extension of the muon profile and the fact that CORSIKA removes particles from the 
longitudinal profile which hit the observation level. 



The air density for a given slant depth decreases with increasing 6. The pions 
decay when their decay length becomes comparable to their interaction length. 
The latter is inversely proportional to the air density and thus pions at the 
end of the hadronic cascade tend to decay at higher energy at larger 6. Muons 
inherit 80% of the energy of their parents on average. Their production energy 
therefore increases from 20 GeV to 100 GeV as 6 increases from 60° to 90°. 

The angular spread of the muons up to this point is mainly caused by the 
transverse momentum pt inherited from the parental pions and the decay an- 
gle between the pion and the muon [7] . An additional small random deflection 
is a caused by the kinematics of the decay. Both effects scale as the inverse 
of the energy. The radial offset of the pions from the shower axis is only of 
the order of a few 10 m and does not contribute significantly to the lateral 
distribution of the muons observed on the ground at r > 100 m [13]. 

After their production, muons are affected by ionization and radiative energy 
losses, decay, multiple scattering, and geomagnetic deflections. Below 100 GeV, 
the energy loss is mainly due to ionization, about 2 MeVg~^cm^ [T9j, which 
translates to about 2 GeV (50 GeV) for a shower at 6* = 60° (90°). The decay 
length is proportional to the energy, for a 10 GeV muon it is 66 km. Due to the 
increase of the production energy with 6, the decay length of the average muon 
is always larger than the distance d from the production point to the ground 
(see Fig. [2]). This explains why a significant fraction of the muons reaches the 
ground at all zenith angles. 

Multiple scattering in the electric field of air nuclei randomizes the directions 
of muons to some degree and erases small scale correlations in the lateral 
profile of the muon density n^. The effect remains small for the total angular 
divergence of the muons from the shower axis up to about 6 ^ 80° where 
multiple scattering becomes the dominant source of the angular divergence 
apart from geomagnetic deflections. 

Geomagnetic deflections 6x for muons in air showers can be approximated 
as [Tllig 

2E^/c 

where e is the elementary charge, d the distance between the muon production 
point and the ground along the shower axis, E"^ is the muon energy, and i?T 
is the perpendicular component of the geomagnetic field with respect to the 
muon direction. Without the geomagnetic field B, the air shower development 
is symmetric in 0. The dependency of the perpendicular component B'^iO^cf)) 
on breaks this symmetry. 

In the extreme 10 GeV muon at a zenith angle of ^ = 60° (80°) gets 

a lateral displacement of about 40 m (1700 m) after a propagation distance 
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Fig. 4. Average divergence of muons from the shower axis at ground level, normalized 
to the average angle between the parent pion and the shower axis, as a function of 
the zenith angle of the shower. The normalized pion-to-muon angle (solid line) is 
determined by the kinematics of the decay and independent of 6; the normalized 
angle due to multiple scattering (dashed line) depends on the accumulated slant 
depth up to ground; the normalized magnetic deflection (dotted line, assuming 
the maximum effect in Malargiie) is proportional to the square of the traveled 
geometrical distance. 

of 10 km (66 km) to the ground level. The impact on the shape of the lateral 
profile is significant. Still, since 5x ^ d, the total number of muons A^^ on the 
ground does not depend on the azimuth angle 0. 

Fig. m shows the relative magnitude of all effects that contribute to the lateral 
spread of the muons observed on the ground. The influence of the geomagnetic 
field becomes important at ~ 70° and dominant at 6' > 80°. 

Muons in the early arriving part of the shower travel shorter distances and 
have smaller inclinations with respect to the ground plane than muons in the 
late arriving part. These effects cause an early-late asymmetry which also 
needs to be taken into account if the lateral profile of the muon density on 
the ground is considered. 

The loss of radial symmetry in the shower front plane due to the geomag- 
netic deformation at very large inclinations and the early-late asymmetry are 
illustrated by Fig.O 
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Fig. 5. Lateral profiles of the muon density in the shower front coordinate system 
of the two showers from Fig.[3l a) ^ ~ 60°, b) ~ 82°. The azimuthal angles of 
the showers are such that the geomagnetic deflections are close to the maximum. 
The radial symmetry around the shower axis is still mostly intact in a), but clearly 
broken by geomagnetic deflections in b). The left-right asymmetry corresponds to 
the early-late asymmetry mentioned in the text. 
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Fig. 6. Same as Fig.O but showing the averaged over the polar angle as a 
function of r^/^ in the shower front coordinate system. Different thresholds ap- 
plied to the muon energy show the same approximate functional behaviour 
Inn^ ~ const. +air^^'^. 

Finally, we point out an empirical approach to describe the lateral profile of 
the muon density on the ground in very inclined air showers. We observe 
that 



dN„ 



K 



n„ 



(X exp(-^afe 



(3) 



dxdy 

describes the output of detailed simulations remarkably well. The coefficients 
ctfc depend on the energy threshold E^^'^'^^ of the detectors, but only weakly on 
the properties of the primary cosmic ray. Even with K = 1, Eq. ([3]) is a good 
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Fig. 7. The lines represent the residuals around fits of Eq. ([3]) with = 3 (this 
work) and Eq. ([4]) (NKG-type) to the lateral profile of the muon density from 
Fig. [6] a) with > 0.1 GeV. Eq. ^ approximates the simulated profile better. Its 
residuals remain around 5 % up to distances of 10 m to the shower axis. 

first approximation. Fig. [6] demonstrates this. 

Eq. ([3]) has a structure very different from the classical formulas of the NKG- 
type that are typically used to describe the lateral profile of the muon den- 



whereas the are another set of coefficients. Formulas of the NKG-type 
diverge at r = 0. 

The choice of a NKG-type formula is generally not motivated by a deeper 
theory. Kamata and Nishimura theoretically derived a formula of this form 
for the lateral profile of a purely electromagnetic shower |5j, but the angular 
divergence in their work is based solely on Coulomb scattering. This is a very 
good approximation for electromagnetic showers, but not for muonic showers, 
as shown earlier. 

Eq. ([3]) always remains finite. It fits better to detailed simulations than Eq. (jlj) 
with the same number of coefficients. An example is shown in Fig. [71 The 
accuracy of the approximation can be improved by increasing K. We note, 
that Eq. ([3]) can also be applied well to the lateral profile of the electron or 
photon density. 



sity ^nm\- 



(4) 
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The form of Eq. ([3]) also allows a convenient way of fitting the coefficients a^. 
We will exploit this and the guaranteed smoothness of the lateral profile due 
to the multiple scattering in Section HI 

2.3 Energy scaling and shower universality 

It was discussed in the previous section that the distributions of the muon 
energy and the angular divergence from shower axis are a function of the air 
density at the altitude h and the slant depth Xatm — -^max and distance d from 
production point to the ground. The amount of geomagnetic defiection also 
depends on the distance d. 

In very inclined showers with 60° < 6 < 90°, Xatm is much larger than Xmax- 
Furthermore, X^ax depends only logarithmically on E and A [21 US]- Thus, 
-^atm — -^max, h-, and d Vary only little with the cosmic ray energy in the range 
10^^ eV <E < lO^o eV (see Fig.E]). 

As a consequence, the muon density profile factorizes in good approxima- 
tion into a normalized profile and the total number of muons A*"^ on the 
ground (see also ref. [7l[5| [T^[TB] ) : 

n^(x, e, 0, E, A) ~ N^{e, E, A) x /^(a;, y, 6, for 6 > 60°. (5) 

The normalized profile depends on the momentum distribution of the muons 
at the production point and the propagation effects to the ground. Both are 
are approximately independent of the energy E and mass A of the cosmic 
ray due to the very slow variation of Xatm ~ A^max, h, and d. This invariance 
property is called shower universality. The total number of muons A^^ also 
does not depend on the azimuth angle 0, as stated earlier. 

If the profile is the one seen by an array of shielded detectors with a 
significant muon energy threshold, should be chosen such that it includes 
the inefficiency of the detector to observe low energy muons. By doing so, A"^ 
can still be identified with the physical number of muons that arrive on the 
ground. 

This factorization approximation was already discussed qualitatively [7t[8|[T3] 
based on simulations with the air shower code AIRES and the hadronic 
interaction models QGSJetOl 1^26] and SIBYLL [27] . Since it is such a cen- 
tral property, we confirm the approximation with the simulation code COR- 
SIKA [^ using the hadronic interaction models QGS Jet-II [23] and FLUKA [21] 

Fig. [8] and Fig. [9] show that the normalized profile varies within 5% between 
10^* eV and 10^° eV. A variation of the cosmic ray mass A between the two 



13 



0.2 



parallel to B-p 



perp. to B-p 



0.1 







-0.1 



-0.2 



p, QGSJet-ll,e-83° 

-10^^°eV 

lO'^-^eV 

' ' ' I ' ' ' ' I ' ' ' ' I ' ' ' ' I ' ' ' ' I ' 

-4 -3 -2-10 1 



10^'° eV 
10^^^ eV 
10'°°eV 



2 3 4 
r / km 



Fig. 8. Variation of the normalized density profile of muons at ground level as a 
function of the cosmic ray energy E for proton showers simulated with the hadronic 
interaction model QGSJet-II. The points at r < (r > 0) show the variation a 
V'-segment around -0 = 0° and ^ = 180° {ip = —90° and ip = 90°) in the shower 
front plane coordinate system. The error bars are dominated by shower-to-shower 
fluctuations in the simulations. To obtain each profile, 10 proton showers from the 
library described in Section [3. 21 were first normalized to one at r = 1000 m and then 
averaged. The geomagnetic field effect in the selected showers is maximal. 

extreme cases of proton and iron nuclei also yields a variation of of 5 %. 

The observed variation of can be compared to the shower-to-shower fluc- 
tuations of the total muon number at ground for cosmic protons which 
are at the level of 13 % (3 %) in simulations of proton (iron) air showers with 
QGSJet-II. The intrinsic shower-to-shower fluctuations of are a natural 
limit to the obtainable energy resolution from a surface array of muon coun- 
ters. The systematic impact of the factorization approximation is negligible in 
comparison if the cosmic ray composition is mostly light. 

The theoretical uncertainties in the normalized profile are estimated by 
comparing |16j showers simulated with the hadronic interaction models QGSJet- 
II and EPOS [28]. They are also at the level of 5 %. 

The total number of muons A*"^ produced around the shower maximum de- 
pends strongly on the energy E and mass A of the cosmic ray. To a lesser 
degree, N'^ also depends on the zenith angle 9 since the charged pions decay 
slightly earlier into muons at larger inclinations, as stated earlier. Since the 
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Fig. 9. The plot shows the same as Fig.[8l but this time the energy is kept at 10^^ eV 
and the cosmic ray mass and the used hadronic interaction model at ultra-high 
energies is varied in the simulation. For each proton (iron) profile, 10 showers (2 
showers) from the library described in Section 13.21 were normalized at r = 1000 m 
and averaged. 

average muon energy increases with 6*, N'^ has to decreases because of energy- 
conservation. Again due to the slow variation of d and Xatm — -^max at large 
inclinations 60° < 9 < 90°, the muon attenuation factor a = N^/N'^ depends 
only on 9 in good approximation. 

The form of the energy dependence of can be calculated with simplified 
Heitler-models of the hadronic cascade [15] and turns out to be a power law. 
We can summarize these findings in another factorization: 

N^{9, E, A) ~ d{9) X K{A) x for 9 > 60°, (6) 

whereas /3 < 1 is a very weak function of the energy E and mass A of the 
cosmic ray and K{A) depends only on the mass A. The energy dependency of 
(3 can be neglected in the range 10^® eV < E < 10^° eV. The mass dependency 
is at the level of 2 % for a change from proton to iron cosmic rays [13l[T6] . The 
function d{9) differs slightly from the true attenuation function a{9) since it 
also parametrizes the decrease of A^"^ with 9. 

The theoretical uncertainties in the parameter /3 and the function d{9) can be 
propagated into iV^. They are at the level of 10 % in the range 60° < 9 < 88°. 
The total theoretical uncertainty of is dominated by the uncertainty of the 
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factor K{A), which depends strongly on the details of hadronic interactions at 
ultra-high energies. The predictions for K{A) of different hadronic interaction 
models differ at the level of about 30% - 50% which is of the same order as the 
difference between proton and iron showers within a single model [T ^fTBlE^ . 

We conclude that the absolute scale of the total number of muons on the 
ground is rather uncertain, but the normalized distribution of the muons 
in very inclined air showers is well-defined. 

If an experiment is mostly interested in the reconstruction of A*"^ from sampled 
muon densities with a model of we expect a theoretical uncertainty smaller 
than 5 % according to the discussed effects. A model of should provide at 
least the same level of precision in comparison with full air shower simulations. 



3 Monte-Carlo simulation of muon ground profiles 

The input for the phenomenological model of the muon density on the 
ground is a library of simulated air showers. 

In this section we shortly present the full Monte-Carlo simulations of very 
inclined air showers that we use in this article to derive the n^-model. Further, 
we demonstrate a fast simulation approach that can be used to speed up 
the simulation of very inclined air showers and in particular generate a large 
amount of input for the model of in a short time. 

3. 1 Theoretical uncertainties and statistical weight- sampling 

The modeling of hadronic particle interactions in extensive air showers is sub- 
ject to considerable theoretical uncertainties. Most hadronic interactions in 
the shower are soft processes with small momentum transfer. So far, such in- 
teractions cannot be calculated within the fundamental theory of quantum 
chromodynamics. 

Calculations of such interactions are based on effective theories and phe- 
nomenology which need to be fitted to data gathered in accelerator exper- 
iments. The center of mass energy in the collision of a 10^° eV proton with 
a nitrogen nucleus is about 400 TeV and far beyond the range of the avail- 
able data. It therefore requires a bold extrapolation of these models. This 
explains the rather large theoretical uncertainties in contemporary air shower 
simulations which showed up mainly in the factor K{A) of Eq. (|6]). 

A technical problem in detailed air shower simulations is the large number of 
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Table 1 

Parameter distribution of the proton shower hbrary generated with 
CORSIKA/QGSJet-II/FLUKA [22H11]. The showers in the Hbrary are dis- 
tributed in smah finite regions of the parameter space. The random distribution 
within each region is given in the table. Each region contains five showers. There 
are 1800 showers in total. 

Parameter distribution low edge of region width of region 

lg£;/eV flat 18.0, 18.5, 19.0, 19.5, 20.0 0.1 

9/° sin(6l) cos(6l) 60,70,74,78,82,86 2 

ct)/° flat 0, 30, 60, 90, 120, 150, 180, 210, 10 

240, 270, 300, 330 



secondary particles. A cosmic ray proton of 10^° eV produces about 10^^ sec- 
ondary particles. Most of them are photons, electrons, and positrons with low 
energies. To calculate such a shower in reasonable time, a statistical weight 
sampling method [5U] is applied during the shower simulation, commonly 
called thinning. This procedure introduces additional artificial fluctuations to 
physical observables of the shower but conserves the mean values. The thin- 
ning phase is usually made as short as possible in order to keep the artificial 
fluctuations small compared to the natural fluctuations that appear during 
the shower development. 



3.2 Full simulation 



The analyses presented in this article are based on a set of 1800 proton showers 
which were generated with the CORSIKA ^2] shower Monte-Carlo code. For 
the high and low energy interactions, the QGSJet-II [23] and FLUKA [2l] 
models were used, respectively. 

We made a choice to distribute the showers in small finite bins in the three- 
dimensional parameter space of energy E, zenith angle 9 and azimuth angle 
0. This leaves gaps where no showers were simulated while other parts get a 
higher statistical coverage. Tabled] summarizes the distribution of the showers. 
Technical details on the simulation options, the used thinning parameters and 
low energy cut-offs can be found in Appendix [X] 

The calculation of one shower takes roughly Iday on a Pentium Xeon 2.8 GHz 
CPU. 
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3.3 Fast simulation 



The geomagnetic field B breaks the 0-symmetry of the simulation of very 
inclined air showers and therefore the detailed simulation has to be run for 
each angle 0. Furthermore, each of these simulation spends a lot of CPU time 
on the calculation of the electromagnetic cascade. 

If only the muon component of the shower is of interest, it is not necessary 
to follow the bulk of the low energy electromagnetic particles. Photons and 
electrons with high enough energies are able to produce further pions in elec- 
tromagnetic interactions with air nuclei. The process with the lowest energy 
threshold of about 0.2 GeV is the production and decay of a A-resonance. 
All electromagnetic particles with lower energies can be dropped from the 
simulation. 

We can further exploit that the development of the hadronic cascade is par- 
tically independent of the geomagnetic field B. The propagation distance of 
the hadronics in the field up to the shower maximum is only a small fraction 
of the propagation distance of the muons and the hadronic cascade develops 
mostly develops at very high energies. Therefore, the geomagnetic deflection 
can be neglected for hadrons. 

This means that it is sufficient to simulate the development of the hadronic 
cascade only once for each combination of energy E and zenith angle 6 up 
to the point where the muons are produced. At this point the simulation can 
be interrupted. The muon positions and momenta can then be rotated to all 
desired values of and the detailed simulations of the muon propagation to 
the ground started from this point. 

This idea was explored with the air shower simulation code AIRES [25J and 
a custom muon propagation code. AIRES was modified in order to write out 
the position, energy, direction, and statistical weight of the muons at their 
respective production points, provided that their probability P = e~'^^^"'^'^'^^ to 
reach the ground is not negligible. High energy thresholds for electromagnetic 
particles make sure that the time-consuming calculation of the electromagnetic 
cascade is stopped early, see Appendix [Bl 

The muons are then rotated to the desired direction of the shower. It is op- 
tionally possible at this point to reduce the weight-sampling of the muons in 
order to increase the statistical precision of the ground profile. Since the cylin- 
drical symmetry of the shower is still intact at the production point, a muon 
with weight w can be replaced by n clones, with each a weight w' = w/n. A 
random rotation around the shower axis is applied to each clone. The cloning 
introduces some artificial correlations to the shower profile at the muon pro- 
duction altitude, but these are erased by random multiple scatterings when 
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the muons reach the ground. The choice of n is an uncritical trade-off between 
the statistical quality of the ground profiles and the calculation time. Values 
between 5 and 50 seem to be reasonable. 

The custom muon propagation code implements the same physical processes 
for muons as the large air shower simulations codes CORSIKA or AIRES. We 
used a custom code to study the muon propagation effects more easily, but 
it is also possible to use either CORSIKA or AIRES to propagate a stack of 
input particles to the ground. 

This method improves the calculation speed dramatically. The calculation of 
the hadronic cascade up to the point where the muons are produced takes 
about 40 minutes on a Pentium Xeon 2.8 GHz CPU. The transformation and 
propagation of the muons for each azimuth then takes only a few minutes 
on the same machine. If 30 values of (p are used for each pair {E, 6), this leads 
to a total speed up factor over 1000 with respect to the simulation of the full 
showers while yielding practically identical results. 



4 Phenomenological model of muon ground profiles 

We have seen in the previous section that it is possible to factorize the ground 
profile of the muon density into a normalized density profile and the 
total muon number iV^ on the ground (see Eq. ([5])), where is approximately 
independent of the cosmic ray energy E and mass A. In this section, we will 
parametrize N^^ and separately as a function of the cosmic ray energy E, 
zenith angle 6, and azimuth angle 0, using a the proton shower library from 
Section 13.21 in order to derive an empirical model of the muon density on 
the ground. 

The parametrization of A*"^ is Eq. 1^. The parametrization of is more com- 
plex and done in two steps. In the first step, the normalized profile of the 
muon density of each individual shower is parametrized. In order to do this 
the normalized muon density is sampled in the shower front plane coordi- 
nate system. The empirical relation between Ig and a/t demonstrated by 
Fig. E] suggests an poljTiomial expansion of Ig in a/t. The effects of early-late 
asymmetries and the geomagnetic deflections resemble dipole and quadrupole 
terms, respectively, which suggest a Fourier expansion in ip. The sources of 
deviation described in Fig. H] sum up to a smooth function. This allows one to 
cut of the expansions at low orders and still get a very good approximation of 

It should be emphasized that this approach can preserve all modulations and 
asymmetries found in the simulation, if enough input data is available. How- 
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10 10^ 

d/km 

Fig. 10. Simulation of the total muon number on the ground A^^ as function of 

the distance d between the shower maximum and the ground along the shower 
axis. Points of different colors are taken from different energy intervals, the average 
energy in each interval is shown at right side of the plot. The continuous lines are 
projections of a two-dimensional fit to the simulation (see text). 



ever, the coefficients do not have a direct physical interpretation: they just 
encode the physical information in a convenient way. 

The expansion performed in the first step leads to a set of Fourier coefficients 
for every simulated shower, with a smooth dependence on the direction of the 
shower (^, 0). It is thus possible in a second step to approximate the evolution 
of each coefficient by an expansion in 9 and 0. This will be detailed in a further 
subsection. 

While we focus on the parametrization of n^, we also want to point out that 
the same scheme may be applied to other local obscrvablcs. For example, if 
the response of the ground detectors depends on the energy of the muons or 
on their incidence angle, these quantities or their relevant combination can be 
parametrized in the same way. 
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4-1 Parametrization of 



The total number of muons on the ground A^^ can be easily extracted from 
each simulated shower. As nuclei give the same profile within a multiplicative 
factor, we regard only cosmic protons [A = 1), so that A*"^ is only a function of 
zenith angle 6 and the energy E of the cosmic ray, as shown in Eq. ([6]). A good 
approximation for the functional form of the energy dependency is already 
given in Eq. ([6]) and thus the only open point is the form of the attenuation 
function a{6). 

A good approximation to a{6) is found empirically by drawing A"^ as a function 
of the distance d{6) between the shower maximum and the point where the 
shower hits the ground in log- log scale, as shown in Fig.[Tni The relation is 
close to a line so that one arrives at the empirical formula: 

Ig N^ = Di + D2 \g{d{e)/km) + (3 lg(^/10i^ eV), (7) 

where Di, D2, and (3 are constants. For the set of proton showers simulated 
with QGSJet-II and FLUKA we get: 

Di = 7.3043 ± 0.0066 
D2 = -0.8240 ±0.0036 
/3 = 0.9352 ±0.0020. 



4-2 Parametrization of the normalized muon profile f^ 



The parametrization of the normalized profile of muon density is now ex- 
plained in technical detail. The procedure starts with a set of simulated show- 
ers and follows four steps. 



A) Coordinate transformation and density calculation. Shower Monte- 
Carlos usually provide weighted muons whose positions are described in ground 
plane coordinates. For each shower, the ground coordinates are projected onto 
the shower plane coordinate system as described in Fig.[TJ To optimize the nu- 
merical precision of a polynomial expansion, it is convenient to choose a vari- 
able ranging symmetrically around zero. Furthermore, we choose u in such a 
way that the relation between Ig and u is almost a straight line. Eq. (|3]) 
suggests: 



'r — \ r ■ 
u = 2^ ^^-^ 1. 

'"max V '"min 



While r ranges between rmin to fmax, u ranges between —1 to 1. 
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The shower plane is then divided into 30 bins in and 30 bins in u. The 
normahzed density in any cell is obtained as a sum over the muons from the 
ground particle file falling in this cell: 



/.Ml) = ^^^ with A,,,, = -^A^{rl,-r^), (9) 

^V^i/lcell ^ COS y 

where rj,rj+i are the edges of the interval in r corresponding to the bin in u 
defined above, and A*"^ is the total (weighted) number of muons in this shower. 

B) Local parametrization. The logarithmic profile Ig/^^ of the shower is 
now parametrized by a polynomial expansion in u and a Fourier expansion in 
ip up to third order: 

3 3 3 

Ig Ur, ij) = J2u'x(Y Ckj cos(j^) + Y Skj sin(j^)) (10) 

fc=0 S=0 j=l ^ 

The parametrization can be fitted to the sampled profile Ig/^ or a shower 
with the linear least-squares method, reducing the fitting problem to a sim- 
ple matrix inversion, and yielding coefficients which are statistically unbiased 
and have minimum variance. At orders higher than 3, the coefficients become 
statistically compatible with zero. 

C) Global parametrization. Each parameter C^j or 5*^^ of the previous 
step is now regarded as a function of (6*, 0) and parametrized using the whole 
sample of simulated showers. The dependence in is a Fourier expansion up 
to the fifth order. The dependence in 6' is a polynomial of the fifth order, again 
conveniently described in the range [d^m-, ^max) through the reduced variable 

V = 2{6 — Omin)/ (^max — ^min) " 1: 



CkAO, 0) = E ^"^ f E C'.^mi ^osm + Y: S',,^, sin(£0) 

m=0 \e=0 1=1 / 

Sk,{e, ct^) = Yv-\Y c'l.^i cos(^0) + E s'l.^i Mm 

m=0 \l=0 1=1 / 



The total number of parameters used to fit the normalized profile is 28 x 
66 = 1848. This number may seem large, but only so because the model 
is 4-dimensional and the number of parameters per dimension multiply. The 
number of free parameters per dimension is small and thus the model is ac- 
tually quite predictive. The parametrization is compared to one of the input 
showers in Fig.lTTl 
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Fig. 11. Example of a) a simulated normalized profile of the muon density in 
shower front plane coordinates and b) the corresponding result produced by the 
parametrized model. The model reproduces the physical structures very well while 
ignoring the statistical fluctuations present in the simulation. 

4.3 Precision of the pammetrization 

The parametrization procedure is applied to the 1800 CORSIKA proton show- 
ers described in Table[Il The bias and an estimate of the precision of the model 
is obtained from an analysis of the residuals of the total number of muons on 
the ground (A^^ — N^)/N^ around the model and the distribution of the resid- 
uals of the normalized lateral profile (/^ — f^)/ f^. 

The analysis of the A^^^-model is straight forward. The bias of the model over 
the whole dataset is negligible. The precision of the A'^^-model is obtained by 
dividing the data into subsets with varying cosmic ray energy and direction 
and analyzing the bias ((A^"^ — N^)/N^) of the model in these subsets. A part 
of the analysis is shown in Fig.[T2]a). A precision better than 2 % is observed 
for the A'^-model. 

For the analysis of the /^-model, the shower front plane is again divided into 
a grid in the coordinates r and ip- The grid is optimised for each individual 
shower in such a way that each cell contains at least 400 explicitly tracked 
Monte-Carlo particles with varying weights. The muons in a local cell carry 
roughly the same weight and thus this procedure limits the statistical fluc- 
tuations o"stat within in the simulated air showers to about 5 %. Cells which 
contain particles with r < 150 m are excluded from the analysis, because these 
particles are additionally thinned in the input simulation with a probablity 
oc [22]. Such cells contain particle with large and strongly varying weights 
so that the amount of statistical fluctuation cannot be estimated in this simple 
fashion. 
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Fig. 12. Shown are a) the average of the residual distribution of the total muon 
number A'^ on the ground as a function of the zenith angle 9 and b) the residual 
distribution (solid line) around the parametrization of the normalized muon density 
in the range 150 m < r < 4000 m. Also shown in b) are the mean and standard 
deviation of the distribution and a Gaussian fit (dashed line) for comparison. 

The residual distribution is shown in Fig. [121 The parametrization shows a 
very small bias of 1 % which is caused by fitting Ig instead of as input 
and the fact that (Igx) 7^ ^g{x) for any random x. The data shows a relative 
fluctuation of 9 % around the model. This observed fluctuation cr(/^) is the 
quadratic sum of the systematic effect cxmodei caused by inaccuracies in the 
model and the statistical fluctuations Ustat in the raw data: 

^'(/.)^^Ldcl(/.)+4at(/.)- (12) 

Substracting (Jstat ~ 5 % from the observed spread yields a precision of the 
/^-model of about 7 %. 

The simulated muon density = A*"^ x is therefore reproduced by the 
model with a global precision better than 2 % and a local precision better 
than 7 %. The latter is only slightly worse than the precision of the shower 
universality approximation from Section 12.31 of about 5 %. To achieve such 
a precision for a variable with the huge dynamic range of about 8 orders of 
magnitude is quite remarkable. 



5 Conclusion 



In this article, we have presented a practical procedure to derive a phenomeno- 
logical model of the ground profile of the muon density from a large set 
of very inclined simulated air showers with zenith angles between 60° and 
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88°. The model is based on a general paramctrization of which exploits 
the smoothness of and the empirical observation that depends on the 
radial distance r from the shower axis approximately like exp(— ay^). The 
parametrization scheme can be adapted to profiles of other ground observ- 
ables. 

As an example, the parametrization procedure was applied to a set of 1800 
proton showers simulated with the hadronic interaction models QGSJet-Il and 
FLUKA in the ultra-high energy range 10^^ eV to 10^° eV. The derived n^u- 
model shows an overall bias less than 2 % in comparison to the simulation 
input. The local precision is better than 7 % and very close to the maximum 
possible accuracy of 5 %. The latter is the level of validity of the shower 
universality assumption in the regarded energy range, which is one of the 
foundations of the model. 

We have further demonstrated a way to speed up the detailed Monte-Carlo 
simulation of the muon component in air showers by a large factor, keeping 
the whole information on the muons at ground level. The speed-up is achieved 
by discarding most of the calculation of the electromagnetic component and 
exploiting the azimuthal symmetry of the muon component around the shower 
axis at the production point. 

A natural application of the n^^-model derived in this article is the recon- 
struction of the energy of very inclined ultra-high energy cosmic rays from 
data collected by a ground array of particle detectors, like the Pierre Auger 
Observatory. The high precision and the flexibility of the model reduce the 
systematic uncertainty of this reconstruction. 
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A Full simulation: Thinning and energy thresholds 



For the correct simulation of inclined showers, CORSIKA ^22j| is compiled with 
the CURVED and UPWARD options. The SLANT option is not mandatory 
for the analyses presented in this article, but necessary if also the longitudinal 
particle profiles of the shower are to be used. 

The showers are thinned according to the strategy described in ref. [3l]. The 
energy E'thin where the thinning starts and the allowed maximum weight i^max 
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for a single Monte-Carlo particle are proportional to the primary energy E of 
the cosmic ray. 




IQ-^'E 
IQ-^E/GeV 
10-^E/GeV . 



The energy dependent Wmax(-E) makes sure, that the amount of actually cal- 
culated particles is roughly independent of the primary energy E, which is 
then also true for the computation time of the showers. 

As a quality/computation time trade-off, a lower maximum weight for hadrons 
and muons w^^^^^ is chosen than for electrons and photons w^^, since muons 
dominate in inclined showers. 

The simulation drops particles below certain momentum thresholds. These 
thresholds have to be adapted to the detection thresholds of the regarded 
detectors. The following cut-offs are used 



for photons, electrons, muons, and hadrons, respectively. 



B Fast Simulation: Thinning and energy thresholds 

To generate the muon distribution at their production point for the fast simula- 
tion explained in Section [373| the following thinning and momentum thresholds 
are used in the air shower code AIRES [25j: 



Pthres 
hadron 
thres 



250 keV/c 
250 keV/c 
0.1 GeV/c 
0.1 GeV/c 




5 X 10"^^ 



0.2 GeV 
0.4 GeV 
0.1 GeV 



(B.l) 



imeson 



0.15 GeV 



E 



'^thres 
• ■ ■ • 

'thres 



0.125 GeV . 



28 



